Stationary, Axisymmetric Neutron Stars 
with Meridional Circulation in General Relativity 
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We present the first stationary, axisymmetric neutron star models with meridional circulation in 
general relativity. For that purpose, we developed GRNS, a new code based on a fixed point iteration. 
We find a two-dimensional set of meridional circulation modes, which differ by the number of vortices 
in the stream lines of the neutron star fluid. For expected maximal meridional circulation velocities 
of about 1000 km/s, the vortices cause surface deformations of about a percent. The deformations 
depend on the shape of the vortices close to the surface and increase with the meridional circulation 
velocity. We also computed models of rotating neutron stars with meridional circulation, where 
neither the surface rotates nor does the rotation velocity exceed the circulation velocity. 
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I. INTRODUCTION 

The shape of a neutron star depends on the motion of 
the stellar fluid. This motion can be decomposed into 
two components with respect to the rotation axis, which 
contains the center of mass and points in the direction 
of the total angular momentum. The first component is 
differential rotation around the rotation axis. This kind 
of motion, which leads to a flattening of the neutron star 
surface, has been investigated intensively, also in general 
relativity [1, 2]. The second component is meridional cir- 
culation, where the flow occurs in meridional planes, i.e. 
planes containing the rotation axis. Its influence, which 
may play a role for new-born neutron stars, is not yet well 
understood. Earlier investigations involve severe limita- 
tions, as they were either restricted to the Newtonian 
framework [3] or they were of perturbative nature only 
[4]. We are interested in how meridional circulations of 
arbitrary strength may influence the shape of a general 
rclativistic neutron star. Moreover, contrary to the ear- 
lier studies considering non-perturbativc fluid motions, 
where in each case only one of the two kinds of motions 
was considered alone, we also investigate the combined 
effect of differential rotation and meridional circulation. 

For that purpose, we developed a new program, called 
GRNS (^Generally Rotating Neutron Star), as the com- 
plexity of Einstein's field equation makes a numerical ap- 
proach to our investigation unavoidable. GRNS combines 
the methods applied in two programs created by two of 
us more than a decade ago. The first one is RNS [1, 2] , a 
general relativistic code for computing models of rapidly 
rotating neutron stars. This code is based on a fixed 
point iteration, and the theoretical considerations of Ko- 
matsu et al. [5]. However, RNS is limited to differentially 
rotating configurations and does not allow for meridional 
circulation. The situation is just opposite for the sec- 
ond program (E. Miiller [3]), which computes meridional 
circulation in non-rotating Newtonian configurations. It 
uses a Newton-Raphson iteration scheme to solve for the 
stationary equilibrium configurations, and describes the 



flow by means of a scalar function, the so-called stream 
function. Below we will show how the stream function 
method can be extended to general relativity. We do not 
use the approach of Komatsu et al. [5] to describe the ef- 
fects of meridional circulation on the curvature of space- 
time, but rely on the more general, theoretical investiga- 
tion of Gourgoulhon et al. [6] , who formulated Einstein's 
field equation as a set of covariant Poisson equations. The 
latter can be rewritten as Poisson equations in flat space 
such that Green functions can be used to solve them. 
We have incorporated these theoretical extensions into 
GRNS, including the fixed point iteration of RNS, and 
a modified Tolman-Oppcnheimer-Volkoff solution as the 
initial guess. That way, GRNS is able to handle differen- 
tially rotating neutron stars with meridional circulation. 

GRNS uses a simplified neutron star model. Consid- 
ering the neutron star during a sufficiently short time 
interval after its formation, it is an appropriate approx- 
imation to assume stationarity. In addition, we limit 
ourselves to axisymmetric configurations. For simplic- 
ity, the stellar matter in our simulations is assumed to 
be a perfect fluid described by a barotropic equation of 
state, although the theoretical formulation presented in 
this paper is kept general enough to deal with arbitrary 
equations of state. We further assume that the neutron 
star has no substructure (i.e., crust, core, etc.), and that 
its chemical composition is homogeneous. Finally, we ne- 
glect the influence of magnetic fields. 

The paper is organized as follows. In Sec. II, we 
rewrite the covariant Poisson equations of Gourgoulhon 
et al. [6] to flat-space Poisson equations and extend the 
Newtonian stream function method [3] to general rela- 
tivity. Intermediate steps during these computations are 
listed in the Appendix. Sec. Ill gives a closer look at the 
GRNS code explaining the fixed point iteration scheme 
and showing how the flat-space Poisson equations can be 
solved by means of Green functions. Actual neutron star 
models are presented in Sec. IV, and we conclude our 
paper in Sec. V. 
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II. THEORY 
A. Notations and conventions 

We use geometrized units for the theoretical compu- 
tations, where the speed of light c and the gravitational 
constant G are set equal to unity. The results obtained 
from neutron star simulations are presented in the more 
appropriate cgs units. The general relativistic formalism 
adopts Einstein's sum convention and is kept as close 
as possible to the conventions chosen in Gourgoulhon et 
al. [6]. However, we make a slightly different choice for 
the range of tensor indices. In our specific choice of co- 
ordinates, concrete indices for 4-tensors, 3-tensors and 
2-tensors are used as follows: 



a,b, 
m, n, 



g {t,r,e,<f>}, 
g {r,e,4>}, 
g {r,e}, 



(i) 



with the time t, the radius r, and the angles 9 and </>. 
Following Misner et al. [7], we choose (—,+,+,+) as the 
signature of the metric. 



Using a special choice of coordinates, which are denoted 
as MTCMA (maximal time slicing - conformally minimal 
azimuthal slicing) in Gourgoulhon et al. [6] , these authors 
write the 2-metric k mn such that the only non- vanishing 
components are k rr = A 2 and kgg = A 2 r 2 , where A shall 
be called the 2-conformal factor. Adopting the definitions 
[6] 



v = IniV, 
a = In A, 



3 = ln> 



M 



r suit 

we then call the eight functions 



v,a 1 f3 1 M r ,M e 1 N r ,N e ,N' 1 ' 



(2) 
(3) 

(4) 



(5) 



basic geometry fields. In contrast, for circular spacetimes 
(i.e. without meridional circulation), only four geometry 
fields are required. 

The neutron star fluid is described by a stress-energy 
tensor 



T a p = (e + p)u a u p + pg aj3 , 



(6) 



B. Basic fields 

We begin our theoretical considerations by defining a 
set of physical fields, the so-called basic fields, which de- 
termine the state of the neutron star uniquely. Due to 
the assumption of stationarity and axisymmetry, these 
and all other fields depend on only two coordinates, the 
radius r and the polar angle 9. There are two types of 
basic fields. The first set describes the curvature of space- 
time, and its fields are denoted basic geometry fields. The 
second set are hydrodynamic fields, called basic matter 
fields, which specify the properties of the stellar fluid. 
In the following, we will define the basic geometry fields 
without giving an explanation of how these fields come 
into being. For that purpose, we refer to Gourgoulhon et 
al. [6]. 

The curvature of spacetime is described by the metric 
g a p. Using a (2 + I) + I decomposition of spacetime, 
the ten independent components of the 4-metric g a $ are 
rewritten to new fields in Gourgoulhon et al. [6]. These 
fields consist of two scalars fields, the 3-lapse N and the 
2-lapsc M , two vector fields, the 3-shift N a and the 2- 



shift M m , and the 2-metric k 



( 9tt 


9tb 1 




V 9at 


9ab J 





, respectively 
N C N C - N 2 

-N a 



-N b 

h a b 



with the 3-mctric 



and 



N a 



-M n M 2 + M P MP 



h ab N h 
h M™ 



with the total energy density e, pressure p, and 4- velocity 
u a . Due to the constraint 



1. 



(7) 



the 4-velocity u a possesses only three independent de- 
grees of freedom. One of them is the specific angular 
momentum, i.e. the ^-component of 



/ _ a 

1M 



(8) 



Under our symmetry assumptions (Sect. HE), the other 
two degrees of freedom can be expressed in terms of a 
stream function ip, an approach known to be working fine 
in Newtonian gravity (equations (2) and (3) in Eriguchi 
et al. [3]). Further below (Sect. HE), we will demon- 
strate that this approach holds in general relativity, too. 
There, we will also explain how the 4-velocity u a can be 
obtained from the stream function ip and the specific an- 
gular momentum 1$. However, by convention we do not 
call the stream function ip a basic field, but the modified 
stream function [3] 



Xo 



r sin# 



So, we refer to 



(9) 



(10) 



as basic matter fields. Actually, we have chosen the 
twelve basic fields (5) and (10) in such a manner that 
they vanish for an empty Minkowski spacetime, which is 
advantageous for our numerical method. 
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C. Geometry equations 

Each of the basic fields (5) and (10) is specified by a 
field equation. In this section, we will address the basic 
geometry fields, which are governed by Einstein's field 
equation 



G a p = 87rT Q f 



(11) 



where G a p is the Einstein tensor. Unfortunately, this 
equation is written in such a compact manner that it is 
practically impossible to directly solve for the basic ge- 
ometry fields (5). Therefore, we use the theoretical con- 
sideration of Gourgoulhon et al. [G] , who rewrite equation 
(11) as a set of eight covariant Poisson equations, each of 
the form 



A cov $ = S, 



(12) 



where A cov , $, and S are a covariant Laplacian, a po- 
tential, and a source, respectively In that paper, the 
covariant Poisson equations are labelled (B3-B7). Try- 
ing to validate these equations, we found that the three 
equations (B3), (B4a) and (B4b) contain minor errors. 
In Appendix A, we show the corrected versions of these 
equations. 

Let us now have a closer look at the covariant Lapla- 
cian A cov appearing in equation (12). This operator is 
the scalar product of two covariant derivatives. Such 
derivatives can be split into a partial derivative and 
corrections arising from the curvature of spacetime, ex- 
pressed in terms of the connection coefficients. In a simi- 
lar manner, it is possible to express a covariant Laplacian 
A cov in terms of the corresponding flat-space Laplacian 
Afl a t and curvature corrections. That way, equation (12) 
can be rewritten as 

A flat $ = S\ 

with some new source S' . In this way, we arrive at the 
following flat-space versions of equations (B3-B7) in [6]: 



3 Ais 

lb 



s v 



2 A[rsin(9(/3 + z/)] 

e 2(a+i/) j^n 

2 A(a + u 



' V ) 

j A a b A^ 6 = Sfr, 



Sp, 

cm 

S a . 



(13) 
(14) 
(15) 

(16) 

(17) 



The flat-space Laplacians 2 A, 3 A, 2 A m „ and 3 A a b are 
listed in Appendix B, and the rather lengthy sources S v , 
S a , Sp, S™j and S% in Appendix I, respectively. Modi- 
fying the sources Sp, S 7 ^, and S a properly, it is possible 
to rewrite equations (15-17) such that the Laplacians in 
these equations act directly on the basic geometry fields 
j3, M m , and a, but not on expressions containing these 
fields. However, the numerical stability of the fixed point 
iteration is highly sensitive to such changes, and in that 
case the iteration would diverge. By experimenting with 



the GRNS code, we found that the choice made in equa- 
tions (15-17) produces convergent results. 

For rotating neutron stars, we achieved convergence 
only by setting the source equal to zero on the grid 
lines adjacent to the rotation axis. Even when increasing 
the grid resolution, only the adjacent grid lines have to 
be modified, and thus setting the source equal to zero is 
a valid procedure. 



D. Matter equations 

The basic matter fields (10) are fixed by the equation 
of state and the equation of general relativistic hydrody- 
namics 



VpT afi = 



(18) 



where V Q is the covariant derivative. As will become ap- 
parent later (Sect. II G), for a stationary configuration, 
the equations of hydrodynamics are integrable under the 
assumption of a barotropic fluid, where the total energy 
density e is a function e (p) of the pressure p only. 

Similar to Einstein's field equation (11), equation (18) 
is again too compact to be solved directly. To overcome 
this problem, we introduce the projector 

= 8% + u a up 

orthogonal to the 4- velocity u a , where Sp is the Kro- 
necker symbol. Then, we project equation (18) in the 
direction parallel and orthogonal to the 4-velocity: 

u a VpT aP = 0, 
qlV[jTa p = Q) 

and bring these two equations in the usual form 



V a [{e+p)u a ] = u a V a p, 
{e + p)u^pu a = ~q afi Vpp, 



(19) 
(20) 



which are the energy equation and the general relativistic 
Euler equation. In the subsequent sections we will extend 
the Newtonian stream function method [3] to general rel- 
ativity by a further reformulation of equations (19) and 
(20). 



E. Stream function 

The stream function tp is introduced to solve the energy 
equation (19). In Appendix C, we show that due to our 
symmetry assumptions equation (19) can be written as 
the vanishing flat-space 3-divergence 



8 i 7 



+ 2 -f 



of the 3- vector 



cot 9f = 



QU", 



(21) 



(22) 
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with [5] 



g = A 2 e 1 (e+p) u t , 
7 = P + v 



(23) 



(note that the symbol g is different than the rest mass 
density p introduced later). Then, similar to the Newto- 
nian stream function method [3], equation (21) is auto- 
matically satisfied by the relation 



1 ( dei' 
gr 2 sind \ -d r ip 



(24) 



This relation reduces two degrees of freedom of the 4- 
velocity u a to the one of the stream function ip (which has 
the units [g/s]). In the Newtonian limit, where u t = —1, 
these two degrees of freedom correspond to the two com- 
ponents u m . Unfortunately, in general relativity, the sit- 
uation is somewhat more complicated, because the com- 
ponents u m together with the third degree of freedom 
also appear on the right hand side of equation (24), hid- 
den in the quantity ut- However, using equation (24), the 
constraint (7), and the definition (8) it is straightforward 
to show that the 4- velocity u a can still be computed from 
the stream function tp and the specific angular momen- 
tum I A,. 



F. Specific angular momentum 

Taking stationarity and axisymmetry into account, a 
short calculation shows that the general relativistic Eulcr 
equation (20) can be written as 

d a p + u a u m d m p 



1 



P 



u /3 u~ f d a gi3 1 — u m d m u a . (25) 



We recall that this equation is a projection orthogonal 
to the 4-velocity u a and has thus only three independent 
degrees of freedom. In the following, we will therefore 
consider only the three spatial components of equation 
(25) by setting a = a in that equation. However, we use 
the temporal component a = t, shown in equation (C3), 
to rewrite these spatial components as 



d a p 
e+p 



1 



u> 3 u' y d a gf3~ t + u t u m d m l a . 



(26) 



Let us start with the azimuthal component of equa- 
tion (26). For that purpose, we set a = (f> and being 
aware that all azimuthal derivatives d^... vanish due to 
axisymmetry, we obtain 



o. 



(27) 



which extends equation (12) of Eriguchi et al. [3] to gen- 
eral relativity. To solve equation (27), we discern the 
three cases 

(A) u m = everywhere, 

(B) u m somewhere, 

(C) u m ^ everywhere. 



In case (A), for which equation (27) is satisfied auto- 
matically, the neutron star can be differentially rotating, 
but without any meridional circulation. We are not in- 
terested in this case, because it has already been inves- 
tigated with the RNS code [1, 2, 8]. Case (B) allows a 
meridional circulation, but not everywhere in the neutron 
star. In this paper, we do not investigate such solutions 
nor do we analyze whether they exist at all. Instead, we 
focus on case (C), where meridional circulation is present 
everywhere in the star. For that case, we use equation 
(24) and rewrite equation (27) as 

deiidrl,!, - d r i>d e l<p = 0. 



This equation is satisfied by 

l<t>=L (VO , 



(28) 



where L (-0) is an arbitrary function of the stream func- 
tion [9]. 



G. Poisson equation for modified stream function 

We proceed with the meridional components of equa- 
tion (26), which can be written as 



— ; \--O m In (u t ) =rwu t 

e + p 2 



utu^dmU (29) 



according to Appendix D, where the quantity if is defined 
in equation (D2). This expression extends equations (7) 
and (8) of Eriguchi et al. [3] to general relativity. Our 
limitation to a barotropic equation of state allows us to 
write the left hand side of equation (29) in the form of 
a gradient. To this end, we introduce the heat function 
[10] 



H(p) 



dp' 



(30) 



For the right hand side, we use equations (24) and (28). 
Then, equation (29) reads 



On 



H{p)+ l -\n(u t f 



w 



gr sin ( 



+ u*L'(V0 )d m i>, (31) 



which implies that the right hand side must be a gradient, 
too, for obtaining a first integral of motion. This leads 
to the condition 



w = gr sin 9 



/WO 



(32) 



with an arbitrary function f (ip). Examining equations 
(D2), (8), and eventually (24), we realize that the quan- 
tity w contains second derivatives of the stream function 
Actually, as shown in the Appendix E, it turns out 
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that there is a Laplacian hidden in that quantity, and 
that condition (32) is equivalent to 



1 



5 A (xo cos < 



(33) 



where the modified stream function xo (9) is used, and 
the source S xo is given in Appendix I. This flat-space 
Poisson equation differs from that presented in Sect. II C 
by the appearance of the angle <f>. 



H. 



Pressure 



We now return to equation (31) and insert condition 
(32). Then, by integration we obtain 



\\n{u t f + H{p)+I^) = C, 



(34) 



with the ancillary function 



i wo 



and the integration constant C. The heat function H (p) 
appearing in equation (34) is invertible, because the total 
energy density e and the pressure p are both positive 
quantities, and thus also the integrand of equation (30), 
such that the heat function H (p) is strictly monotonous. 
Denoting the inverse heat function as (p), equation 
(34) the pressure can be expressed as 



max 



M, 




■max l 



FIG. 1: Initial guess of basic field xo- Each one of the 
four panels shows the distribution of the field xo inside of 
the neutron star for one choice of the pair (M r , Me) accord- 
ing to equation (38). The top, left panel visualizes the case 
{Mr, Me) = (!,!)• Proceeding to the right increases the 
value of the parameter Me, and we have to go down to raise 
the value of M r - For each panel, the maximal absolute field 
value is called max. The values max and —max are repre- 
sented by the brightest red and green colors, respectively. 



that the field xo becomes non-zero. For that purpose, we 
consider surface-adapted coordinates 



R(ey 



(37) 



H( Pc ) 



In 



Ui 



with the central pressure p c 



(35) 

the central stream function 



ipc, and the central covariant temporal component u£ of 
the 4-velocity, respectively. In general, the above inver- 
sion is performed numerically. 



III. NUMERICS 

A. Fixed point iteration 

In Sect. II, we have introduced the basic fields 

v, a, p, M r , M 9 , N r ,N 9 ,N\p, e, X o, h (36) 

and their governing equations (13-17), (28), (33), (35) 
and the equation of state. We solve these equations 
by means of a fixed-point iteration after having speci- 
fied an initial guess vq, ■■■,l ( pfi, which is constructed in 
the following way. We set the basic fields (36) to the 
Tolman-Oppenheimer-Volkoff solution. This spherically 
symmetric solution has no meridional circulation, and 
therefore the basic matter field xo vanishes. Then, we 
modify the Tolman-Oppcnhcimcr-Volkoff solution such 



where R (9) is the radial coordinate of the neutron star 
surface. Then, the invariance of scalars tells us 



and we set 
Xo (r,t 



Xo[f,0) = xo ( 



Xo" aX sin (Mrfn) sin (M g 9 



(38) 



with an arbitrary constant Xo lax an d parameters 
M r , Me e {1, 2, ...}, as shown in Fig. 1. 

To describe the iteration steps, we denote the values 
of the basic fields (36) after the s-th iteration step by 
v s , i^ jS , and those after the previous s-th iteration step 
by v s -i, l<p tS -i- The basic geometry fields are then 
iterated via the equations 



= A~ S u 



2 A 1 S a s -i 



(39) 
(40) 

& = ( 2 A~ 1 S' ( 3 iS _i) / (rsin#) - v a , (41) 

(43) 



A Q i 



jv.s-n 



where Sy a _i denotes the sources of equations (13-17) 
computed from the basic field values v s -ii •••) l<j>,s-i- The 
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inverse Laplacians •••A" appearing above are deter- 
mined with Green functions according to the detailed 
description given further below (Sect. IIIB). Subse- 
quently we evaluate the source of equation (33) using 
the newly computed values v s , . . . , MJ of the basic geom- 
etry fields and the old values p s -i, l<p. s -i of the basic 
matter fields, the outcome being called S XOtS -i- Next, 
we calculate 

Ps = H- 1 [H ( Pc ) + (in K^x) 2 - In K s _i) 2 ) /2 
+J(^ C ) - I(ips-i)} , 

e s = e(p s ) , 

Xo = ( 3 A _1 (5 Xo , s _icos0)) /cos0, 
l$ >a = L (ips-i) , 



(44) 



where s _ 1 is the value of the quantity u£ computed 
from the values of the basic fields v s -\, l<f>,s— x- 



B. Green functions 



In the following, we address the inverse Laplacians 
2 A _1 , 3 A" 1 , 2 A m n 1 and 3 A a b 1 appearing in equa- 
tions (39-44). Each of these inverse Laplacians comes 
from an equation having one of the forms 



2 A$ = S, 
3 A$ = S, 
2 A™,$' 1 = S n , 
S a , 



3 ao 



A Q b $ b 



(45) 
(46) 
(47) 
(48) 



where and S " are a potential and a source, respec- 
tively. To solve for the potential of equation (45), (46), 
we have to compute the first, second, ... of the fol- 
lowing integrals, shown in the same order as above [5]: 



$(£) = 

$ m (r,6>) = 

$ a (r,6>) = 
where 



d x 



S(x>) 



2n/ln\x-x'\ 
dV s & ■ 



47T | X — X 1 1 



(49) 
(50) 



d (r, 9) m 
d (x, z) n 

d(r,e, 



d{x 1 y,z) 



,d(x',z') n S°(r , ,g / ) 
C d(r',0')° 2n/\n\x- x'l 1 

3 , d(x , ,y,'z') b S c (r',6') 
X a(r',6»',0') C 47r|x-x'| ; 



x = (x, z) = (rsin#, 7-COS0) (51) 

for the 2-dimensional and 

x = (x, y, z) = (r sin 9 cos r sin 9 sin 0, r cos 0) 

for the 3-dimcnsional integrals (analogous for x') together 
with the Jacobian determinants d (...) a / d (. ..) . The dots 



appearing in equation (50) represent additional terms re- 
quired for the boundary condition of the basic matter 
field xo an d are given in Appendix F. In that Appendix, 
the main focus is on an expansion of the above integrals 
in terms of trigonometric functions and Legendre poly- 
nomials. Typically, only a few terms of the expansions 
arc needed, and neglecting the remaining ones leads to a 
strong reduction of the computational cost necessary to 
evaluate the integrals. 



C. Slicing conditions 

According to Gourgoulhon et al. [6], the basic geome- 
try fields M m and iV a satisfy the following slicing condi- 
tions 



(N 2 M m ) 
N a 



where the 2-and 3-covariant derivatives '||' and '|' are 
defined in Appendix B. A lengthy, but elementary cal- 
culation shows that these two slicing conditions can be 
rewritten as the vanishing flat-space 2- and 3-dimensional 
derivatives 

1 



of the quantities 



d m M t 
2 



-Ml 



0. 



Ml 



n: + cot on: 



e 2(a+u) j^m 



N a . 



Contrary to an analytic one, a numeric evaluation of the 
Green functions of Sect. IIIB always produces values 
for the fields M m and N a which somewhat violate the 
slicing conditions. Therefore, we consider a 2- and 3- 
dimensional Hclmholtz decomposition and set the gradi- 
ent parts of the two quantities M™ and N£ equal to zero 
such that the divergences of the remaining parts vanish 
as demanded by the slicing conditions. Such a procedure 
is valid, because when the fixed point iteration has con- 
verged, the gradient parts (which have to be set to zero) 
decrease with increasing numerical accuracy, i.e. with 
increasing grid resolution and number of trigonometric 
functions and Legendre polynomials used in the integral 
expansions of the Green functions. 



D. Final gauge 

Even now, the basic fields are not yet determined com- 
pletely. There are still remaining gauge degrees of free- 
dom left, because we can add an arbitrary constant to 
each of the potentials appearing in equations (45-48) 
without violating these equations: 



const'". 



(52) 
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For the 3-dimcnsional Poisson equations, the Green func- 
tion 



vanishes in the limit \x — x'\ — > oo. Thus, we choose the 
potentials of these Poisson equations to vanish at infinity 
by setting the constant appearing in equation (52) equal 
to zero. However, for the 2-dimcnsional Poisson equa- 
tions such a choice is not possible, as the Green function 

2 G(x,x' ) = —]n\x-x'\ 

is not bounded for \x — x'\ — > oo. Therefore, we proceed 
as follows: 

The potential r sin (/? + v ) of the 2-dimensional Pois- 
son equation (41) obeys the Dirichlct boundary condition 
according to Appendix F. Hence, we set the constant in 
equation (52) equal to zero in this case, which completely 
fixes the field (3. In case of the 2-dimensional Poisson 
equation (40), we choose the constant in equation (52) in 
such a manner that 

a (r = 0) = (r = 0) . 

For the remaining 2-dimcnsional Poisson equation (42), 
we have to consider only the Cartesian z-componcnt M| , 
because the ^-component obeys the Dirichlct bound- 
ary condition, and we impose the condition 

Ml (r = 0) = 0. 

E. Circulation modes 

Eventually, we have a closer look at the function / (ip) 
introduced in equation (32). We follow the same ap- 
proach as Eriguchi et al. [3] and limit ourselves to the 
power law 

f(i>) = krP n , (53) 
with some constant k and the exponent 

n = 0, 1. (54) 

For the case 

/ WO = ty, 

we find different meridional circulation modes i}) m (with 
m = 0, 1, 2, ... ), as in the Newtonian case [3]. The dis- 
tributions of the field xo belonging to these modes look 
similar to the ones displayed in Fig. 1 but are somewhat 
deformed. 

Unfortunately, when nearly having reached one of the 
higher circulation modes ipi,i/)2i— during the fixed point 
iteration, the fixed point iteration method always starts 
to converge to the fundamental mode tpo- In order to 



obtain the higher modes, we therefore projected the lower 
ones away. For that purpose, we assume that we have 
already evaluated the first m c — 1 modes, i.e. we know ipm 
for m = 0,1,..., m c — 1. Then, the fixed point iteration 
leads to the m e -th mode by replacing 

m e — 1 

lj>-*1p- ^ C m 1pm 
m— 1 

at every iteration step, with adequately chosen coeffi- 
cients C m . If an orthogonality relation 

/ dr deVhW m (r, 9) ip m <ip m > = S mm * (55) 
Jo Jo 

with the weight function W m (r, 8) exists, the coefficients 
C m are given by 

J °°dr f*d9VhW m (r,d) Wm. 
io°° dr So d^VTiWm (r, 0) ij m ip m ' 

where h = det/i a b. However, we neither know whether an 
orthogonality relation exists nor do we know the weight 
functions W m (r, 9) . After some experimenting, we found 
that the choice 

W m (r,e) = e (56) 

is sufficient to achieve a convergence to higher modes. 
This does not necessarily mean that (56) is the correct 
weight function, but it could be very close to it. 

In addition to the usage of (56), we perform the follow- 
ing steps in the GRNS code: The pressure distribution of 
the solutions investigated in this work is always equato- 
rially symmetric. However, in our treatment, equatorial 
symmetry is not guaranteed exactly due to the limited 
numerical accuracy. Therefore, equatorial asymmetry 
may increase during the fixed point iteration, eventually 
leading to an undesired meridional circulation mode. In 
order to prevent this, we symmetrize the pressure dis- 
tribution at every iteration step. A similar method is 
performed for the basic field %o > which has either even or 
odd parity depending on the considered mode. 

F. GRNS 

We have implemented GRNS under Linux in C++, and 
it possesses an OpcnGL visualization interface, which al- 
lows the user to supervise the fixed point iteration. The 
user has full control over the iteration, which can be 
stopped and restarted at any time. The user can se- 
lect any of the physical fields either when the iteration is 
stopped or even when it is running, and display it on the 
screen. It is also possible to visualize the flow of the neu- 
tron star fluid in real time and to see how the neutron 
star surface changes during the iteration. These code 
features are of advantage when analyzing the stability of 
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GRNS, and they were also very helpful in debugging the 
code. 

We used a numerical grid of 150 radial and 156 angular 
zones to compute the models discussed in Sect. IV. For 
the sums arising from the expansion of the integrals (see 
Appendix F), we chose the upper limit 



*=... i=... 



(57) 



where l max = 10 for all models of Sect. IV. This limits 
the computational cost, and we are able to compute in- 
dividual circulation modes within about a minute on a 
current single core CPU. In Appendix G, we discuss the 
convergence behavior and consistency tests of GRNS. 



IV. RESULTS 

A. Reference model 

The results presented in the following were obtained 
with the usual polytropic equation of state 



P 
e 



Kp l 



1 



where K is the polytropic constant, p the rest mass den- 
sity, and r the polytropic exponent, respectively. From 
this equation of state its is readily seen that 



P 



r - i 



(58) 



We keep the maximum absolute value of the stream func- 
tion tp fixed at the value '(/'max such that it does not mat- 
ter which value we choose for the constant k in equation 
(53). The fundamental mode of each neutron star model 
is then unambiguously specified by the following param- 
eters 

max ' 

i.e. by the exponent appearing in equation (53), the 
central pressure and the central density (which fixes the 
polytropic constant K) of the initial guess, the polytropic 
exponent, the maximum absolute value of the stream 
function, and by the distribution (28) of the specific an- 
gular momentum. As a reference model, we choose a 
non-rotating neutron star with parameters 

n = 1, 

p c = 9.1 • 10 34 erg/cm 3 , 
Pc = 7.9-10 14 g/cm 3 , 

r = 2, 

</>max = 3-10 34 g/ S , 

L(i>) = 0, 



circulating at its fundamental mode. This set of param- 
eters corresponds to a neutron star with the following 
properties: 



R = 15.6 km, 
M = 1.51 M , 
Scire = 1043 km/s. 

The average radius R is the average proper radius 



(59) 



R c = / dre a 
Jc 

averaged over all angles 9 (over all radial paths C). We 
compute the rest mass via [11] 

M = / drdOd^VhN^p, 
Jv 

where V is the volume of the star and the average circu- 
lation velocity is given by 



(v r f + (rv 6 )' 



The basic fields of the reference model are shown in 
Fig. 2. The 3-lapse logarithm v is one of the two con- 
tributions (the other one comes from the 3-shift N a ) to 
the time dilation, which is strongest at the center of the 
neutron star. The 2-conformal factor logarithm a repre- 
sents the stretching of space in meridional planes and is 
also strongest in the center. We do not display the basic 
geometry field j3 in Fig. 2, because its values are not sig- 
nificantly different from those of the field a. The 3-shift 
N a represents the dragging of spacetime caused by the 
flow of the neutron star fluid. Since the reference model 
is circulating, this leads to the vortex visible in the lower 
left panel of Fig. 2. The total energy density e is also 
not shown in that figure, because it is correlated to the 
pressure p by the trivial analytic relation (58). Note that 
in the lower right panel of Fig. 2 the outermost contour 
is the surface of the neutron star. As the values on the 
axes of Fig. 2 do not refer to proper but to coordinate 
distances, the average neutron star radius suggested by 
the contour is somewhat smaller than the actually correct 
value (59), which takes the curvature of spacetime into 
account. The 2-shift M m and specific angular momen- 
tum vanish for the reference model, because it does 
not rotate. Instead of considering the modified stream 
function xo = (rsin0), we show the stream function 
tjj itself in the top left panel of Fig. 3, as the contours 
of the stream function ip are stream lines of the merid- 
ional flow. Note that the stream function of the reference 
model contains only a single vortex. 

The distributions of the 3-lapse u, the 2-conformal fac- 
tor a and the pressure p given in Fig. 2 do not show 
any significant difference when compared to the Tolman- 
Oppcnhcimcr-Volkoff solution used to construct the ini- 
tial guess for the fixed point iteration. The only relevant 
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FIG. 2: Basic fields of reference model. The four panels 
show the 3-lapse logarithm v, the 2-conformal factor loga- 
rithm a, the 3-shift TV" and the pressure for our reference 
model. The color coding of the three scalar plots is the same 
one as in Fig. 1, where the maximum absolute field value max 
is displayed at the top of each panel. For the vector plot, the 
quantity max denotes the maximum vector length (note that 
all vectors lie within the paper plane). 
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FIG. 3: Circulation modes - Part 1. The panels show the 
stream lines of the first four meridional circulation modes of 
the reference model. The color coding is the same as in Fig. 1. 
Note that red and green refer to clockwise and counterclock- 
wise motion, respectively. The kinks visible at the neutron 
star surface are a result of the finite numerical resolution and 
the use of contour plots. 



B. Higher circulation modes 



changes are that the 3-shift N a and the stream function 
tjj, which vanish for the Tolman-Oppenheimer-Volkoff so- 
lution, display the distributions shown in the lower, left 
panel of Fig. 2 and the upper left panel of Fig. 3, re- 
spectively. 



The other three panels in Fig. 3 show higher circula- 
tion modes with additional vortices. Actually, each vor- 
tex in the stream function ip is accompanied by a vortex 
in the 3-shift N a (not shown for the higher modes). Let 
us now compare the four modes shown in Fig. 3 with 
the four initial guesses of Fig. 1. With this in mind, it 
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does not matter that both figures show different fields, 
because the difference between the stream function ip and 
the modified stream function \o is merely a factor r sin 9. 
Such a factor causes deformations of the vortices, but 
their number and orientations remain unaffected. Then, 
it seems as if the top left initial guess of Fig. 1 leads to 
the top left mode of Fig. 3, and so on. However, there 
is no such one-to-one correlation. Instead, it can, for in- 
stance, occur that the four modes of Fig. 3 are produced 
by the four initial guesses M r = 1 and Mg = 1,2,3,4. 
Actually, it does not matter too much how the initial 
guesses are chosen. It is only important that they are 
somehow different, and to improve the fixed point iter- 
ation, we have made a choice that is at least similar to 
the expected outcomes. 

In total, we are able to compute 16 modes for the ref- 
erence model (including the fundamental mode) before 
GRNS fails due to too large numerical errors, i.e. we find 
a much larger set of modes than in Eriguchi et al. [3] . Fig. 
4 shows some additional higher modes. The lower right 
panel of that figure represents the highest mode obtained 
where numerical errors do not yet have a significant im- 
pact on the shape of the vortices. Looking at Figs. 3 and 
4, it is obvious that the circulation modes constitute a 
2-dimcnsional mode set. 



C. Surface deformations 

A closer look at the surface of the neutron star in Fig. 
5 reveals that surface deformations of about half a per- 
cent of the average neutron star radius are present. The 
shape of these deformations depends on the vortices in 
the vicinity of the surface. Vortices deeper inside of the 
neutron star have only a small impact on the deforma- 
tions. In Eriguchi et al. [3], the numerical resolution was 
too small to resolve such properties. However, it was 
determined that for the fundamental mode with only a 
single vortex the meridional circulation causes the neu- 
tron star to become slightly prolate. 

The strength of the deformations depends also on the 
velocity of the fluid in the vortices, and thus on the max- 
imum absolute value ^i max of the stream function. To 
test the limits, we increased the value of tpmax in GRNS 
as far as possible. The choice ^ max = 3 • 10 34 g/s for 
the reference model is already close to what is possible 
with GRNS, and we were only able to increase that value 
to about V'max = 10 35 g/s, above which numerical errors 
grow dramatically. However, even at the limiting value, 
the surface deformations do not exceed about one percent 
of the average neutron star radius. 

D. Circulation velocities 

For the choice ipmax = 10 35 g/s, the average circula- 
tion velocity of the reference model is u c i rc = 3431 km/s, 
whereas for a very slowly circulating model with "0max = 




02468 10 02468 10 12 
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FIG. 4: Circulation modes - Part 2 

3 • 10 34 g/s the average circulation velocity is only w c i rc = 
10 km/s. Actually, it is a good approximation to say that 
the average circulation velocity v c scales linearly with the 
maximum absolute value ?/> max of the stream function. 
However, for the average radius and the mass of the neu- 
tron star, we do not find any significant changes when the 
circulation velocity changes, within the range considered 
here. 

Tab. I lists the average velocities v c i rc of the circula- 
tion modes. For that purpose, we use a mode numbering 
similar to Fig. 1: the first mode number in Tab. I gives 
the number of vortices in x-direction, and the second one 
in z-direction. The velocity v clIC obviously increases with 
the mode numbers, which is a consequence of the orthog- 
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FIG. 5: Surface deformations. The figure shows the radial 
coordinate R (9) of the neutron star surface depending on the 
angle 9. From thick to thin, the three lines refer to the upper 
left and upper right mode of Fig. 3, and the upper left mode 
of Fig. 4. 

TABLE I: Average circulation velocities. The table 
shows the average velocities v c h c for the eight circulation 
modes of Figs. 3 and 4. 

Mode (1,1) (1,2) (1,3) (1,4) (2,1) (2,2) (2,3) (3,3)' 
v ciTC [km/s] 1032 1296 1569 1778 1364 1661 2002 2444 

onality relation (55). As the average circulation velocity 
Wcirc rises with the maximum absolute value ip max of the 
stream function, changing the orthogonality relation (55) 
by inserting an appropriate factor D m in the integrand, 
it is possible to create new modes with the same velocity 
v c i rc , but a different maximum ij) max . Because it is not 
easy to determine the factor D mi our method is more 
practical. 

E. Constant / (tp) 

Eriguchi et al. [3] investigated both cases (54). There- 
fore, we now take the reference model and set the param- 
eter n = such that we obtain 

/ W = k. 

For this case, we do not find a collection of circulation 
modes but only a single solution, which is shown in Fig. 
6. There, we see that compared to the upper left panel of 
Fig. 3 the stream function is slightly deformed. In addi- 
tion to that, we find that the average circulation velocity 
becomes 

Vdrc = 1366 km/s, 

while the average proper radius R and the mass Ai of 
the neutron star are the same as for the reference model. 
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FIG. 6: Constant f (ip). Stream function tp for the case 
n = 0. The color coding is the same one as in Fig. 1. 

F. Rotating neutron star 

In contrast to Eriguchi et al. [3], we have also com- 
puted rotating neutron star models with meridional cir- 
culation. For that purpose, we set the specific angular 
momentum L (-0) to a non- vanishing value. The simplest 
choice would be 

L (tp) = const ^ 0, (60) 

where the neutron star rotates rapidly near the rotation 
axis. Eriguchi et al. [3] exclude such models, because 
the rotation velocity becomes infinite when approaching 
the rotation axis in Newtonian configurations. In gen- 
eral relativity, the situation is more complicated, because 
spacetime itself can be dragged by the neutron star fluid. 
We are not able to present stringent analytical reasons 
which exclude condition (60), but simulations obeying 
that condition do not seem to converge for higher nu- 
merical resolutions and have to be considered as invalid 
neutron star models. 

Therefore, we consider the next simplest case 

L (tp) = const • ip. 

For the choice 

LW=3 .10 14 -^ — , (61) 

we find the average circulation and rotation velocities 

Scire = 1168 km/s, 

v rot = 431 km/s, (62) 
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where the latter is the average value of 



Vrot 



= r sin 0ip 



When compared to the reference model, the average ra- 
dius and the mass do not show significant changes, just 
as most basic fields. Two exceptions are the specific an- 
gular momentum 1$ and the 2-shift M m , which vanish for 
the reference model, whereas for the rotating model (61) 
their distributions look like as displayed in Fig. 7. They 
now exhibit a vortex, similar to the stream function ip 
and the 3-shift N a . In contrast to the reference model, 
the 3-shift N a vector does not only lie within the merid- 
ional plane as shown in the lower left panel of Fig. 2, but 
it has also contributions perpendicular to that plane for 
the rotating model (not displayed in this paper). 

Having investigated model (61), where the circulation 
and rotation velocities are of the same order of magni- 
tude, the question arises what happens when the rota- 
tional velocity is much smaller or much larger than the 
circulation velocity. Since the reference model itself is 
non-rotating, we continuously increased its rotation ve- 
locity from zero to the value (62) by changing the con- 
stant in relation (61). Currently, GRNS fails to compute 
models where the average rotation velocity is larger than 
about the value (62). Actually, one would expect that 
there are at least a few models where the rotation veloc- 
ity is much larger than the circulation velocity, and pos- 
sible improvements in the numerical method could allow 
their computation. 



G. Higher modes with circulation and rotation 

There are not only higher circulation modes for non- 
rotating configurations, like the reference model, but also 
for rotating models. In this case, the difference between 
the fundamental mode and the higher modes is similar 
to what has been discussed in Sect. IV B. All of the 
basic fields containing a single vortex for the fundamen- 
tal mode, i.e. the four fields ip,N a ,l ( j 3 ,M m , exhibit the 
same number of vortices for a certain higher mode. These 
vortices are always roughly located at the same spatial 
position for the four fields. Similar to the non-rotating 
configurations, the modes where both circulation and ro- 
tation are present turn out to be a 2-dimensional mode 
collection. 

Fig. 8 shows radial profiles of the specific angular mo- 
mentum (per unit rest mass) 



3 = 



P 



-U0 



(different from the specific angular momentum 1$, which 
is defined per unit energy) for the fundamental and one 
higher mode belonging to the rotating model (61). There, 
we do not only recognize combined co- and counterrota- 
tion, but also that in parts of the star dj/dr < 0, which 
implies that the flow does not satisfy Solberg's criterion 
for local stability (see [12] for a proof in GR). 
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FIG. 7: Rotating neutron star. Specific angular momen- 
tum 1$ and the 2-shift M m for the rotating model (61). The 
color coding of the contour plot is the same one as in Fig. 1. 
For the vector plot, the quantity max denotes the maximum 
vector length. 



V. CONCLUSIONS 

We computed the first stationary, axisymmetric neu- 
tron star models with meridional circulation in the frame- 
work of general relativity. For that purpose, we have 
developed GRNS, a new code that uses a fixed point it- 
eration method starting from a Tolman-Oppcnheimer- 
Volkoff like initial configuration, similar to RNS, written 
by N. Stergioulas. 

We found meridional circulation modes like Eriguchi et 
al. [3] in the Newtonian framework. However, by using a 
fixed point iteration instead of a Newton-Raphson one, 
we were able to automatize the process of computing such 
modes. As a result, we identified a much larger set of 
modes than these authors. Our study shows that the 
circulation models form a two-dimensional set. 

The circulation modes differ by a varying number of 
vortices in the stream lines of the fluid. These vortices 
cause deformations of the surface of the neutron star. 
The deformations are influenced strongest by the vor- 
tices in the vicinity of the surface, and their influence 
rises with the value of the circulation velocity. We found 
surface deformations of the order of one percent for ex- 
pected maximal circulation velocities of about 1000 km/s. 
However, the radius and mass of the neutron star do not 
significantly depend on the circulation velocities. 

In contrast to Eriguchi et al. [3] , we also computed ro- 
tating neutron star models with meridional circulation. 
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FIG. 8: Specific angular momentum. The two panels 
show radial profiles of the specific angular momentum j for 
two different modes of the rotating model (61) at angles 8 — 
(thinnest curve), 9 = j, and 6 = \ (thickest curve), respec- 
tively. The upper panel belongs to the fundamental mode, 
displayed in Fig. 7, and the lower one to one of the higher 
modes. The specific angular momentum 1$ and the stream 
function rj) of this higher mode exhibit a vortex distribution 
as in the upper, left panel of Fig. 4. 



For such models, the rotation velocity is highest in the 
center of the vortices, and vanishes at the surface of the 
neutron star. In addition to that, we were unable to 
find rotating neutron star models, where the rotation ve- 
locity is significantly larger than the circulation velocity. 
We are not sure whether this is possibly caused by our 
symmetry assumptions. 

There are clear perspectives for a future application 
of the outcomes of this investigation. Perturbing the 
obtained modes, a dynamical evolution of the neutron 
star could show the influence of meridional circulations 
on gravitational waves, for which a direct experimental 
detection is expected in the near future. Another appli- 
cation is investigating the influence of meridional circu- 
lations on neutron star oscillations. Both methods offer 
a way to determine by means of observations whether 
meridional circulations are present in (young) neutron 
stars. At the current stage we are unable to decide how 
common such circulations are in nature, because with 



our study wc did not evaluate stability criteria for the 
circulation modes. 

In the near future, the two most important extensions 
of this investigation will be changing the topology to a 
toroidal one and to include magnetic fields. We have 
already presented some first thoughts in that direction 
in the Appendix finding that the field equation for the 
specific angular momentum is strongly affected by the 
presence of a magnetic field. 
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Appendix A: Corrected geometry equations 

In our investigation, we use equations (B3-B7) of Gour- 
goulhon et al. [6] , which contain minor errors in the three 
equations (B3), (B4a), and (B4b). We found these er- 
rors by thoroughly verifying all the mathematical equa- 
tions in that paper. For equation (B3), this was done by 
hand, and similar to Gourgoulhon et al. [G] we used the 
computer algebra program Mathematica to validate the 
rather lengthy equations (B4a) and (B4b). In the fol- 
lowing, we will give the three corrected equations. Due 
to their length, we will not list them completely here, 
but show only where the corrections appear. Several 
new mathematical quantities are introduced during the 
derivation of the equations, which will not be defined 
here. Instead, we refer to Gourgoulhon et al. [6], whose 
notation and conventions are adopted by us, except for 
the differently chosen indices (1). 

Equation (B3) of Gourgoulhon et al. [6] has the form 

... =4n(E + S:) + K ab K ab + ^-, 

where the term 2m r m e v^ r V Q is missing on the left hand 
side, i.e. the correct equation (B3) reads 



... + 2m r m e v r V = 4tt (E + S») + K ab K ab + —. 

In the last line of the left hand side of equation (B4a) of 
Gourgoulhon et al. [6], there appears the expression 



M 



M r 



Al 



Here, the presence of the 2-lapse M in the squared 
bracket is an error, i.e. the correct equation (B4a) has 
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the form 



and the two less familiar vector Laplacians 



^ [M\ e {» >e - 4a o) - M\ ee ] ... . 

The second but last line of the left hand side of equation 
(B4b) of Gourgoulhon et al. [G] reads 



m 
'~M 



M 7 



A, 



M 



M 



A 



The presence of the 2-lapse M in the squared bracket is 
an error, and the third term —2A t gg/A has to be deleted 
such that the correct equation (B4b) reads 



2 Am _ fm 2 , 



3 A a b = <5 h a3 A 




(Bl) 



with the matrix element 77133 = j.d r +2^r-de. Note that 
equation (Bl) holds only in case of axisymmetry, where 
the azimuthal derivatives of the Laplacians 3 A a b and 
3 A vanish. 
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Appendix C: Energy equation 



Appendix B: Flat-space Laplacians 

In the following, we consider 2- and 3-dimcnsional 
Laplacians. For that purpose, we use the inverse 2- and 
3-metrices k mn and h ab [6], and the Christoffel symbols 
of the second 

2-prn _ imp ti 
1 no ~ " ± pnoi 

1 be — n 1 dbc 

and of the first kind 

r Q ^j 7 = - {dpg ai + d^gpa — d a gp 7 ) . 

These Christoffel symbols allow us to compute the 2- and 
3-covariant derivatives [6] of a tensor T^---\ 

rpm... o rpm... _i_2 pn rpq... . 2pcy rpm... 

1 n...\\p ~ u P 1 n... ' L qp 1 7n... ' ■■■ L np 1 q... '"J 

rpa... o rpa... 1 3 j^a rpd... i 3t^(/ rpa... 

1 b...\c ~ °c 1 b... + l dc 1 b... ' '" 1 6c J d... — ■•• • 

Then, limiting ourselves to a scalar $ and a vector $ a , 
we construct the following Laplacians: 





= 2 A$, 




= 3 A$, 




— 2 A m $™ 


^ C 4>V 


= 3 A a fc $ b . 



For flat-space, we choose 

k mn = diag (l,r 2 ) , 
h ab = diag (l,r 2 ,r 2 sin 2 6) , 

such that a straightforward calculation gives the two well 
known scalar Laplacians 

2 A = d 2 + \d 2 + -d r , 

3 a n.9 1 o9 2 „ cot „ 1 

3 A = d 2 r + -d 2 e + -d r + —d e + , 2 d 2 , 



In this appendix, we show that for our symmetry as- 
sumptions it is possible to write the energy equation (19) 
as the vanishing flat-space 3-divergence (21) of the 3- 
vector (22). For that purpose, we expand equation (19) 
as 

d a [(e + p) u a ] + (e + p) T p pa u a = u a d a p, (CI) 

where = g aS Tsp^ are the common 4-dimensional 
Christoffel symbols of the second kind. Because of [13] 

with the determinant g = dct g a p , and our assumptions 
of stationarity and axisymmetry all temporal and az- 
imuthal derivatives <9 t ... and c?^... vanish. Hence, equa- 
tion (CI) becomes 

d m [(e + p) u m ] + (e+p) u m d m ln y/=g = u m d mP . (C2) 

To compute the determinant g, we use its relation to the 
determinant k = dctk mn (see equation (2.27) of Gour- 
goulhon et al. [6]), the fact that in MTCMA coordinates 
(Sect. II B) k = A 4 r 2 , and the definitions (2) and (3): 

sf^g~ = NAlVk = e 2Q+ V sin (9. 

For the right hand side of equation (C2), we apply the 
temporal component 

u m d m P = -{e + P ) u m d m In in (C3) 

of the general relativistic Euler equation (25), obtained 
by setting a = t in that equation and using stationarity. 
Hence, equation (C2) becomes 

d m [e 2a+1 (e+p) u t u m ] 

+e 2a+7 (e + p) ut {^u r + cot Ou 6 ^ = 

which is equal to equation (21). 
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Appendix D: Meridional components 

Let us consider the general relativistic Euler equation 
(20) rewritten in the form of equation (26). This ap- 
pendix deals with the meridional components of that 
equation, which means setting a = m in equation (26) 
such that 

dmP 1 a ci nc, 1 mil 

— — = -u u^d m g a p + u t u d n l rn . (Dl) 
For the second term on the right hand side, we write 



u n d n l m = rw I I + u n d m l n , 



with 



w = - (dgl r - d r lg) , 



(D2) 



begin the quantity 'to' defined in Eriguchi et al. [3] in the 
Newtonian limit. Then, using the constraint (7) a short 
calculation shows that 



1 



u a u p d m g a p - u t u n d„ 



-d m In Ut + UtU^dm — . 



Hence, equation (Dl) can be brought in the form (29), 
which has a form similar to equations (7) and (8) of 
Eriguchi et al. [3]. 



Appendix E: Special Laplacian 

In the following, we show that there is a Laplacian 
hidden behind the quantity w defined in equation (D2). 
For that purpose, we extend equations (15) and (17) of 
Eriguchi et al. [3] to general relativity. 

The first step is to define the quantity 



;dgip 



D = dfy + 



which is equal to the left hand side of equation (15) of 
Eriguchi et al. [3]. Then, equation (24) allows us to write 

D = sin 9 [-d r (r 2 gu e ) + d e (Qu r )] . (El) 

Next, we use u m = g ma u a , and the decompositions 
(2.9b) and (2.25b) of Gourgoulhon et al. [6] to obtain 
after a short computation 



The only non-vanishing components of the 2-metric k mn 
are k rr = 1/A 2 and k ee = 1/ (A 2 r 2 ) . Therefore, equation 
(E2) gives 

d r (r 2 gu e ) = d r c g - d r (gig) , (E3) 
d g {gu r ) = d 8 c r - d e (gl r ) . (E4) 
By inverting equation (E2) to 

I --(^c -k 



(E5) 
(E6) 



_ urn rvrnnU 

Ut \ Q 



and using equation (24) we also find 

detp 



4 1 

Q 

h = - [ ce 

Q 



r z sm 

d r ip 
sin 9 



Now we combine equations (D2) and (E3-E6) such that 
equation (El) becomes 



D = sin 9 (dgc r — d r ce) — gr sin 9w 



(d r 4> + eg sin 9) d r 



dgip 



c r sin 9 ) d [ 



(E7) 
\ng, 



which extends equation (15) of Eriguchi et al. [3] to gen- 
eral relativity. 

Next, we introduce the quantity [3] 

ip COS (j) 
^ r sin 9 
to write equation (E7) in the form 
r sin 9 ( cos 



D 



r 



(dgC r — d r Cg) — gw COS (j) 



X H 1 COS ) d r 

r r 



+ ~2 (®8X + X cot ^ — rc r COS <j>) dg 



In g 



In the Newtonian limit, the expression in curly brackets 
becomes the right hand side of equation (17) of Eriguchi 
et al. [3], i.e. equation (33) extends equation (17) of 
Eriguchi et al. [3] to general relativity under condition 
(32). 



Appendix F: Integral expansion 

Below, we expand the four integrals appearing in Sect. 
IIIB in terms of trigonometric functions and Legcndrc 
polynomials. 



u m = u t W 



(E2) 



2-scalar 



with 



gu t 
A 2 ' 



M n 



(M p u p + U(j> ) 



N n N a N r, 



N 2 a N 2 



u t 



For the integral (49) we use the expansion 

1 min' (r, r') 
I max' (r, r') 



ln|x — x'\ — In max (r, r') — 



1=1 



(cos {19) cos (W) + sin (19) sin {19')) 
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of Komatsu et al. [5] . Applying the von Neumann bound- 
ary condition 



S(r,n + 9) = S (r, ir 
a short computation shows that 



$ (r, 9) 



dr'r' In max (r,r') / d6'S(r',6') 



-WlcosiW) I dr'r' 111111 £> r > (Fl) 
d9' cos (W')S(r',9'). 



, i min ; (r, r') 



We use this expansion for the basic geometry field a, i.e. 
equation (40). In that case, the potential is $ = a + v 
and the source S — S a . 

For the basic geometry field j3, i.e. equation (41), the 
integral (49) has to be evaluated. The potential is <fr = 
r sin# (/3 + i>), which has to vanish on the rotation axis. 
Therefore, we apply the Dirichlet boundary condition 



S(r,ir + 6) = -5(r,7r- 



and find 



^ Y~[ ' Jo max 1 (r, ?' ) 



, , min' (r, r') 



d(9' sin (W')S(r', 9'), 



with the source S = S 



2. 3-scalar 



For the integral (50) we use the expansion [5] 



(/ — m)\ min' (r, r') 

la; — x'\ 1 — ' e — ' (I + m)\ max i+1 (r, r') 

i=o m=-l y ' K ' 



EE 



■P[" {cos6') P{ n (cos 6) e™^-*'), 

where P™ are the associated Legendre polynomials. In 
case of an axisymmetric source 

S(x) = S (r, 9) , 

it is easy to show that 



$ M) = 4E P '( 



1=0 



dr'r' 2 



min' (r, r') 



max' +1 (r, r') 

f dO'Pi (cos 6') sine'S (r',6'). 
Jo 



We apply this result to equation (39), and choose the 
potential $ = v together with the source 5* = S v . 



The integral (50) also appears in the equation for the 
basic matter field \o (44). In that case, the source is no 
longer axisymmetric but obeys 

S(x) = S(r,0)cos0. 

Then, it is possible to show that 



with 



$o (r,9) 



$ (r, 9, (/>) = $(r, 0)cos( 
$M) = $oM) + - 



1 00 1 

-J^TT, ^P, 1 (COS 9) 

2 ^ 1(1+1) 1 v ; 
i=i y ' 



(F3) 



(F4) 



dr r 



>>2 



min 1 (r, r') 



max 1+1 (r, r') 

[ dO'Pl (cos9')sm9'S(r',9'). 
Jo 



The dots appearing in equation (F3) serve the same pur- 
pose as in equation (50), i.e. they represent a remaining 
degree of freedom for the choice of the boundary condi- 
tion of the field xo- Looking at equation (19) of Eriguchi 
et al. [3] we actually see that 



$ (r, 9) = $q (r, 0) + ^T air 1 Pi (cast 



(F5) 



i=i 



with arbitrary coefficients a;. We have to choose these 
coefficients in such a manner that there is no flow across 
the surface of the neutron star. For that purpose, we 
consider surface-adapted coordinates (r, 6) (Sect. Ill A), 
where the radial coordinate of the neutron star's surface 
becomes f = 1 (37). The condition for no flow across the 
surface is 

u r = 

(the letter 'S' denotes that this relation holds only on 
the surface). Using equation (24), it is straightforward 
to show that this leads to 

d § ^ I 

As 



ip = r sin 9xo = fR ( 9 ) sin 9x 



we then find 



R 9 



tan %o + tan 9d §X o = (F6) 



Though not being the only mathematical solution, we 
currently limit ourselves to the case 



s n 
Xo = 



(F7) 
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which is the Dirichlct boundary condition and which 
leads to 

d§Xo = 

because in surface-adapted coordinates the surface cor- 
responds to a constant radial coordinate f. This shows 
that the Dirichlet boundary condition (F7) implies the 
constraint (F6) for no flow across the surface. In terms 
of the potential 

$ (r,(?) = $(r,0) 
equation (F7) implies 

i (1,0) = 0. 

For this boundary condition, it can be shown that 
If 2/ + 1 pl 



(F8) 



dO'Pl (cos0') sin0'$ o (1,0') , 



with 



$ r,0 = $ o (r,0). 



Note that we use the expansions (F4) and (F8) both for 
the potential $ (r, 0) — \q (r, 0) and the source S (r, 0) = 
S xo (r,0). 



4. 3-vector 

The fourth integral of Sect. MB is used to compute 
the potential $ a = N a and the source S a = S$f- In this 
case, the derivation of the expansion in terms of Legen- 
dre polynomials is somewhat lengthy, but still straight- 
forward such that we give only the result: 



$ a (r, I 



smtf 

cos 6 



Pl (cos0) 



' drV 2 nW (r,r') 



max' +1 (r, r') 

/ d0'F ; 1 (cos0')sin0' 
Jo 



sin0'S r (r',0') +r' cos0'S" e (r',0') 
r' sin 0'^ (r',0') 

/ COS0 



J2 p i ( cos( 

J 



V 



' drV 2 min ( r ' r ') 



max' +1 (r, r') 

/ d0'P ; (cos 0') sin 0' 
Jo 

(cos0'5 r (r', 0') - r' sm9'S e (r', 0')) . 



Appendix G: Tests 



3. 2-vector 

Next, we consider the third integral of Sect. IIIB to 
evaluate the potential $ m = e 2 ( a+ ^M m and the source 
S rn = S™j. Similar to Sect. Fl, we have to specify 
boundary conditions. We use the 2-dimensional Carte- 
sian coordinates introduced in equation (51), in which the 
potential and the source have the components ($ x ,3 ,z ) 
and (S x , S z ), respectively. We apply the Dirichlet bound- 
ary condition on the ^-component <f> x , and the von Neu- 
mann one on the z-component $ z . These two quantities 
are then governed by equations resulting from replacing 

$(r,0) -> ® x (r,0), 
S (r',0') -> S x (r',9') 

in equation (F2) and 

$(r,0) $ 2 (r,0) , 
5 (r',0') S*(r',0') 

in equation (Fl). 



Here, we present convergence and consistency tests 
performed with GRNS. For that purpose, we introduce 
the convergence indicator 

„ Sgrid We — -Fs-l| 

C s = lOOmax— ^— — , (Gl) 

F Z^grid^sl 

where F s is the distribution of the basic field F at the it- 
eration step s, and t The sums extend over the numerical 
grid. First the fraction is evaluated for all basic fields F, 
and then the maximum value determines the convergence 
indicator C s . That way, this quantity is most sensitive 
to the basic field which converges least. Note that the 
convergence indicator becomes C s = for perfect con- 
vergence. Fig. 9 shows convergence tests for three differ- 
ent resolutions. We see that for the fundamental mode 
GRNS converges for all three resolutions, i.e. the con- 
vergence indicator C s approaches zero during the fixed 
point iteration. For higher modes, the convergence indi- 
cator drops initially but then starts to fluctuate, never 
getting close to the value zero. Improving the weight 
function given in equation (56) might improve this be- 
havior. We have also considered different numbers of 
terms Z ma x used in the sums (57), namely 3, 10, 50, which 
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Iteration step s 



FIG. 9: Convergence tests. Dependence of the conver- 
gence indicator C s defined in equation (Gl) on the iteration 
step s for a model similar to the reference model. The three 
solid lines refer to grid resolutions of 59 x 52, 150 x 156, and 
501 x 507 zones (the higher the resolution, the thinner the 
corresponding line). The blue dashed lines refer to the first 
higher mode. 



where F a p = d a Ap — dpA a is the electromagnetic field 
strength expressed in terms of the electromagnetic 4- 
vector potential A a . Moreover, we do not only have to 
solve Einstein's field equation (11), but also Maxwell's 
field equation 

where p q is the charge density. A lengthy, but straight- 
forward computation shows that the general relativistic 
Euler equation (20) becomes 

(e + p) u' 3 W p u a = -q af3 W pP - p q F a p u^ 

and equation (27) changes to 

ma , PqU m (IfdmAt + d m A^j 



results is a convergence behavior similar to that shown 
in Fig. 9, when the three grid resolutions are replaced 
with the three values of ? max . Eventually, we have val- 
idated the correct behavior of GRNS, when the spatial 
extension of the numerical grid is changed. 

Appendix H: Outlook on the magnetic field 

If a magnetic field is included, the stress-energy tensor 
(6) has the additional contribution 



Hence, the simple solution method (28) does no longer 
work if a magnetic field is present. 

Appendix I: Sources 

In Sect. IIC, we have written the covariant Poisson 
equations (B3-B7) of Gourgoulhon et al. [6] to flat-space 
Poisson equations. The sources S v , S a , Sp, SJ}, and 
of these equations, and the source S Xo of the Poisson 
equation (33) are listed in the following. Similar to Ap- 
pendix A, several new mathematical quantities appear 
below, which are defined in Gourgoulhon et al. [6]. 

The scalar sources are 



S v = A 2 J 4tt (E + S a a )+ K ab K ab + ^ - 



and 



and 



r\ 2 



-42 + K") 



(rA) 



+ K) 2 



(i/ e ) 2 - (m r ) 2 v, rr - (m 9 ) i/ ( 



—2m r m 9 v tr g — (m r rn r r + m 9 m r e ) v T — (rn r m 6 r + m 9 m 9 e ) v^g — 2m r ra 9 v^ r v^ — /? r ^ r — 
S Q = A 2 jsvrs + 1 [(q r + U>m r ) K, r + (q 8 + com 9 ) K,g] + \_K r [M, q] r + Kg [M, q} 9 ] + 3n m ^ 



(if l? 17171 I \ T T T 

~ 77 \"-mn"' ~r~ ~r~ J - / mn lj 



{V,r) ~ ~ 
\ r 



Sp = — {SvrMiVs™ - 2«v [M, q] r - 2n [M, q} 6 - M (q r + um r ) n. r - M (q 9 + urn 6 ) k,q 



+MN (K mn K mn + k 2 - L mn L mn ) } - r sin 6 
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and 



1 



S Xo = — (dgc r — d r cg) — ggr sin 8 



The two components of the 2-vector source are the r-component 

\2 at2 



i^r ~ U * L ' w ) + \{ drXQ + v + ^) dr + ^ {deXa + cot 9x0 ~ rCr) 



d e In g. 



S r M = A Z N 2 U{a + v) r M r (a + v) r + M\ r + — {a + v) g M r (a + v) e + M\ 



—W {a + v) e + S" M + 2M r S c 



and the ^-component 



S% = A 2 N 2 \ 4 (a + v) r \M 9 (a + v) p + M 6 ' 1 + -j (a + v) e hi 9 (a + u) 8 + A/ e 



+ ^Af e (a + i/) jr + ^Af (a + v) fi + S"m + 2M e S a J 
The three components of the 3- vector source S 1 ^ arc the r-component 
Sfr = A 2 ^-16nNJ r -2K rr N, r -2K re N,e-2m r m e N r rg ~N r r ^2 



a. r + m r m e {2a. g — 



M 



M 



N r A2 



2(m e ) 2 



(rAY 



— [M\ r + 2A/ 6 



+ 2m r m e {pLfi - a.e) - 2m r 



2N" 



■J2 + ( m ) M-m M 



1 \ m r M r ( 

a.g +m , m ( /x r + 4a, r + - J ' 



A/ 7 ' 



M M 



e . m I 2, 



I 2A/ r M e 



1 (im r ) 
+ - + 

r r 



+ 2mW (M% - M\ r ) - -3 



-l-^) 2 + (r^) 2 



A'/ r 





' 1 


-9- 


A 2 



A/, 



A/ 



+ m r m e (— /i,rA*,i 



4a. <?// r - 2 2a. r a. e + — — 2—— + 2 (m ) 



a. r ) — 

r A 



+ W [2A/F r (/v - 4a, r 



M e r + 2a e ) + 2A/% ( /j r - a r - - ) - 2M r rr - M <rl 



r 
M r 



— [M\ r (^o - 4a, e ) 



r0\ 



M 2 



( \ 2 J. M 



2 (M\ r Y + M r r M e g + r 2 (M r Y\ \ ~ N ° \ jr a ^,r + m r m s ±a fii i fi - 2 (ag) 
A 



M A 
V 



+ 2 f m ) a r a a + 2u e a r — + — 2A4 (/i a - as ) 

rr AIM' 



2M\ g [ 3a, r + - ) + M\ e {ji fi - Aa fi ) - M\ BB - 2M\ r6 



^ [M\ e (p,e ~ 4a, e ) - M r >e0 ] 



w {M\ e M e fi + 2M r r M r g + r 2 M e r M%) \ - (rr/) 2 ^ N\ rr + ( - + M,r ) AT, r 



^ + Ov) 2 



AT' 



the ^-component 



r 



N r -\p fi N\ e + (-l3, g + 2cot6f3, r + (3, rg + 2l3, r p, e ) N°, 



S% = A 2 I -16nNJ e - 2K er N, r - 2K ee N fi - 2m r m e N B rf) - N r r { -2 



(rA) 



+ (m e ) 2 



a g + 2m r m (fx r — a r ) 
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M 



m r _1 _ M 6 



2m & M\ 



-N r , 



a, r + 2(m e ) -2m 6 



(rA) 



A I 



+ 2 (m r ) 2 J a r + m r m e (fx >e + 4a fi ) + — [m r (2M r r + M%) - m 9 r 2 M e r ] \ - N 6 \ e \ 2 



1 



i 

A? 



i 6 ) 2 ] a fi + m r m e ( 2a, r - fi >r + + jj (m r M e >r + m e AI\ r ) J - N* T J (",« ~ /V) 



1 

',4 2 



1 (Aro r ) : 



+ 2m r m e (Af r . r - M %) + 



1 



+ (m r ) 2 — (rm 6 ) 
1 



.}-**{■ 



A 2 



(m r ) 2 + (rm e )' 



1 

A 2 



(m ) n, r fJ,,6 + {m ) -jj- 



(rA)- 



(2cot0/3, r + /3, re + 2/3, r /3, e ) 



(rA)' 



A 



A 



+ 2(m e f 



a H j a rM e + m m 

(rA) 



,/V a 



1 



2^1 _ 4_£ _ + 4« rM _ 2 ( a , r r - (/v) + 



M, 
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A 



rl) 



A 



m 
~M 



fi r — 4a r I M „ — M _ 
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m 
~M 



M 



4a, r - - ) - 6M 9 r a, e - 2M 6 9 ( -/x, r + o, r + - ) - M r >rr , - 2M" >r( ; 



M r r M e r + 2M e >r M%) } - iVV 2 



2a, r a.0 H rr 2- 



+ 2(m e )' 



OLfififi + m r m 



2 2 a r u g + 4a r u 

r r ' 



AI 
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1 \ 2 1 ^,66 


TTl r 






-2(m B f 


(«,«) + ^ 
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2AF 
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M r r M\ 
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(m 6 ) 2 (N\ gg + M , e 7V%) + (m e ) 2 (cot 2 9 + 1- f} M ) N e ) 
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N r 
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